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COLLAPSE  OF  PARTIALLY  MIXED  REGIONS  IN  STRATIFIED  FLUIDS 


Introduction 


There  is  a  sizeable  literature  devoted  to  the  collapse  of 
regions  of  mixed  fluid  immersed  in  density  stratified  media.  The 
physical  application  is  the  description  of  energy  exchange  between 
the  stratified  fluid  and  the  region  of  mixed  fluid.  The  mixed 
region  is  generated  by  an  overturning  internal  wave  (Long  [1])  or 
by  the  turbulence  in  the  wake  of  a  moving  obstacle  (Schooley  and 
Stewart  [2],  Schooley  [3]).  Thus,  the  generation  of  the  mixed 
region  is  a  turbulent  process  and  the  mixing  is  accomplished  by 
the  turbulent  eddies.  The  mixed  region  grows  in  size  by  entrainment 
of  surrounding  fluid  until  the  kinetic  energy  in  the  turbulence  is 
dissipated.  The  mixed  region  is  statically  unstable  and,  when  the 
turbulence  decays  sufficiently,  it  proceeds  to  collapse  to  its 
level  of  equilibrium. 


There  have  been  numerous  attempts  to  model  this  collapse 
process.  Except  for  models  put  forward  by  Ko  [4]  and  Bergin  [5], 
the  collapse  stage  has  been  completely  decoupled  from  the  mixed 
region  growth  stage.  Most  physical  models  have  taken  as  an  initial 
condition  a  two-dimensional  circular  region  of  quiescent  mixed 
fluid.  The  physical  models  have  been  represented  mathematically 
in  numerous  ways  to  describe  the  ensuing  collapse  process.  Schooley 
and  Stewart  [2]  considered  the  mixed  region  to  be  a  small  perturba¬ 
tion  from  the  ambient  equilibrium  state  and  solved  the  linearized, 
inviscid  equations  governing  the  resulting  fluid  motions  by  an 
eigenfunction  expansion.  The  mixed  region  was  assumed  to  be  circular 
in  sh.  De  with  a  density  perturbation  of  Gaussian  form.  This  method 
has  been  applied  again  by  Schooley  and  Hughes  [6]  but  the  mixed 
region  was  assumed  to  be  fully  mixed  in  this  case.  In  either  case, 
the  linearized  governing  equations  are  not  strictly  valid  unless 
the  density  perturbation  is  indeed  a  very  small  one.  The  case  in 
which  the  amount  of  mixing  is  small  but  uniform  across  the  mixed 
region  can  be  validly  treated  in  this  manner  and  it  has  been  solved 
by  Hartman  and  Lewis  [7.. 


The  fully  mixed  case  Is  not  as  interesting  physically  as  the 
original  diffusely  mixed  case  modelled  by  Schooley  and  Stewart  [2] 
because  of  the  sharp  discontinuity  in  density  at  the  edges  of  the 
mixed  region.  On  the  other  hand,  this  discontinuity  is  helpful  in 
distinguishing  the  fluid  interior  to  the  mixed  region  from  that 
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exterior  to  the  region.  The  experiment  of  Wu  [8]  took  advantage  of 
this  discontinuity  by  tagging  the  mixed  region  fluid  with  dye.  This 
idealization  has  been  helpful  also  in  various  mathematical  models 
of  the  process.  Mei  [9],  Padmanabhan  et  al  [10],  and  Bell  and  Dugan 
[11]  have  solved  models  that  account  for  nonlinearity  in  the  govern¬ 
ing  equations  but  which  omit  realistic  interactions  of  the  mixed 
region  with  its  surroundings.  Wessel  [12]  and,  more  recently,  Dugan 
et  al  [13]  and  Young  and  Hirt  [14]  have  solved  the  full  nonlinear 
equations  numerically  for  the  case  of  sharp  discontinuity  in  the 
density  at  the  edge  of  the  mix^d  region.  The  numerical  method  is 
the  most  complete  and  accurate  method  if  proper  care  is  taken  to 
verify  “he  accuracy  of  the  calculations.  Thus,  it  appears  to  be 
the  most  useful  method  for  any  further  studies  of  mixed  region 
collapse. 

The  studies  referred  to  above  have  brought  to  light  most  of 
the  physical  mechanisms  present  in  the  collapse  phenomenon.  This 
understanding  is  quite  limited  though  in  that  knowledge  is  sparse 
for  the  more  realistic  case  in  which  '.he  turbulence  does  not  com¬ 
pletely  homogenize  the  fluid  in  the  mixed  region.  In  particular, 
what  are  the  mixed  region  shapes  and  energy  balances  when  the 
mixing  is  uniform  but  not  complete  and  when  the  mixing  is  not 
uniform.  These  questions  are  taken  up  here  and  the  numerical 
method  utilized  by  Dugan  et  al  [13],  and  modified  as  indicated  in 
the  appendix,  is  helpful  in  answering  them. 

Analysis 

The  mixed  region  is  taken  to  be  circular  in  shape  and  to  be 
at  rest  initially.  The  density  structure  is  linearly  increasing 
with  depth  outside  the  region  so  the  Brunt-Vaisa la  frequency  there 
is  constant  to  the  accuracy  of  the  Boussinesq  approximation.  In 
the  mixed  region,  the  density  profile  is  taken  to  have  several  forms 
as  shown  in  Figure  1. 
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Figure  1 


The  case  I  is  assumed  to  be  uniformly  mixed  in  that  the  density 
gradient  in  the  interior  is  constant  but  is  less  than  that  outside 
the  region.  That  is,  the  density  profile  is 


1-P*  r>CXo 

\  i-M  ri-Q0  (1) 


where  r  is  the  radial  distance  from  the  center  of  the  mixed  region. 
The  case  II  has  the  density  profile 


(2) 


and  the  case  III  has 

F  M  (3) 

The  latter  case  in  a  sense  is  intermediate  between  the  case  I  having 
a  sharp  edge  and  the  case  II  having  a  very  diffuse  edge.  In  the 
nonuni formly  mixed  cases  II  and  III,  the  density  perturbation  is 
continuous  and  the  mixing  is  complete  only  at  the  mixed  region  center. 
The  amount  of  potential  energy  stored  in  the  mixed  region  varies 
with  the  type  of  mixing  and  it  is 

i,  PE  =  \?oNil0-e^ao4 
n.  PE 
III.  PE 


In  each  case,  the  potential  energy  is  defined  as  the  wcrk  required 
to  construct  the  mixed  region  from  the  linearly  stratified  surround¬ 
ings  in  a  reversible  manner.  The  cases  II  and  III  have  the  same  a- 
mount  of  potential  energy  and  that  amount  is  one  half  that  of  a  fully 
mixed  region  of  the  same  ’radius’. 


The  governing  equations  are  the  Navier-Stokes  equations  and, 
in  the  absence  of  turbulence,  the  nonlinear  effects  are  important 
but  the  viscous  and  diffusive  ones  are  not  for  reasonable  scales. 
Several  finite  difference  methods  of  solving  the  nonlinear  equations 
are  available.  Wessel  r  12]  initially  solved  a  problem  of  this  type 
in  which  the  mixing  was  complete  with  a  scheme  similar  to  the  one 
of  Dugan  et  al  r 13  .  The  solutions  obtained  were  in  good  agreement 
with  the  experiment  of  Wu  [8,  but  the  energy  conservation  properties 
of  the  code  were  uncertain.  Young  and  Hirt  [  1A]  solved  the  same 
problem  but  the  method  used  appears  to  dissipate  the  available 
energy  quite  rapidly.  The  modified  MAC  code  used  therein  is  stable 
only  because  of  the  presence  of  nonzero  viscosity  and  the  temporal 
history  of  the  mixed  region  width  does  not  compare  well  with  the 
experiment.  The  method  used  by  Dugan  et  al  [ 13]  has  been  shown  to 
conserve  energy  quite  well  so -it  is  a  useful  one  to  study  the 
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energetics  in  the  problem  at  hand.  It  was  also  shown  to  predict 
results  very  similar  to  those  of  the  experiment  of  Wu  [S]  lor  the 
fully  mixed  case.  The  numerical  method  is  a  modified  form  of  that 
originally  put  forward  by  Williams  [15]  ;  the  formulation  is  laid 
out  in  Dugan  ct  al  [13]  and  it  ha9  been  modified  as  indicated  in  the 
appendix  to  allow  grid  stretching  and  to  better  conserve  energy. 

Uniform  Mixing 

The  case  I  of  partial  but  uniform  mixing  is  considered  first. 

The  numerical  solutions  for  the  mixed  region  shape  versus  time  are  shown 
in  Figure  2  for  £  =  .81  and  .49  where 

e  =  P  ‘  <5> 

and  ^  and  |3  are  defined  in  expression  (1)..  The  solution  for  €  *  .49 
shows  little  change  from  the  fully  mixed  case  in  the  early  stage  of 
motion.  However,  eventually,  all  the  shapes  for  even  those  values 
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of  £  approaching  zero  are  different  from  the  fully  mixed  case.  The 
equilibrium  positions  of  all  fluid  particles  in  the  partially  mixed 
region  are  above  or  below  the  mixed  region  centerline  so  that  the 
shape  never  collapses  to  a  line  on  the  axis  of  the  cylinder  as  it 
does  in  the  fully  mixed  case.  The  final  equilibrium  shape  does  not 
appear  to  be  far  from  an  ellipse  in  cases  for  which  £>,5". 

The  width  of  the  mixed  region  as  a  function  of  time  is  shown 
in  Figure  3  for  several  values  of  £  . 


In  all  cases,  the  width  increases  to  a  maximum  and  exhibits  a 
slightly  underdamped  oscillation  to  its  equilibrium  width. 

This  case  of  uniform,  partial  mixing  is  simple  enough  to  be 
amenable  to  analysis.  Hartman  and  Lewis  [7]  have  obtained  an 
analytical  solution  for  this  case  and,  to  the  validity  of  their 
solution,  the  width  of  the  mixed  region  is  easily  shown  to  be 

al0/a  =  2-e-za-e:)  ^(N-OiN-n1  <6> 
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This  solution  is  valid  only  when  the  amount  of  mixing  is  small;  that 
is,  it  is  valid  in  the  limit  of £- 1  .  However,  it  is  plotted  in 
Figure  4  and  it  compares  well  with  the  numerical  results  for  short 
time  for  all  cases  in  which  computations  were  made.  The  nonlinear 


effects  are  clearly  apparent  in  the  mixed  region  widths  in  Figure  4 
in  the  cases  in  which  the  mixing  is  greater  than  207.  (c4.fi  )  as  the 
numerical  and  the  linear  solutions  digress.  Even  in  the  case  of  €=  ‘SS’ 
though,  the  solutions  are  not  identical.  The  nonlinear  solutions  ' 
outdistance  the  linear  ones  and  they  do  not  oscillate  as  quickly. 


> 


The  solution  method  of  Bell  and  Dugan  [llj  in  which  the  energy 
in  the  mixed  region  is  assumed  to  be  conserved  can  be  extended  to 
apply  to  the  case  of  partial  mixing.  The  governing  equation  that 
corresponds  to  equation  (12)  in  Bell  and  Dugan  [ 11]  is 

(l  +  o<4)  *  Z-  (l-2e^0 

The  solution  of  this  equation  for  small  time  goes  like 


o.&Y, 

/i n 


1  +  ~  nY-  0(N4*4) 


A* 


Of  course  this  expression  is  never  strictly  valid  since  the  mixed 
region  never  conserves  its  energy.  On  the  other  hand,  it  is  not 
restricted  to  €.  -  1  .  Dugan  et  al  [13]  have  modified  this  analysis 
with  the  assumption  that  there  is  an  equipartition  of  kinetic  energy 
inside  and  outside  the  mixed  region  and  have  verified  the  resulting 
analytical  prediction  with  the  numerical  solutions  in  the  case  of 
tomplete  mixing.  Making  that  assumption  here,  the  governing  equation 
is  t 

2(i>of4)i  +(i-2eo0  (1-2.0=  O 

This  cannot  be  solved  in  closed  form  as  in  the  fully  mixed  case  but 
the  solution  valid  for  small  time  is 


1  -  £ 
8 


mV-  o(mV) 


(7) 


It  is  of  interest  to  note  that  the  short  time  behavior  of  the  solution 
of  Hartman  and  Lewis  [?,  (see  equation  (6)  above)  is  identical  to  this 
expression.  The  comparison  between  the  numerical  results  and  the 
asymptotic  solution  (7)  above  is  shown  in  Figure  3. 


The  energy  content  of  the  mixed  region  can  be  obtained 
approximately  from  the  foregoing  results.  Bell  and  Dugan  [ll]  were 
able  to  compute  energies  under  the  assumption  that  the  mixed  region 
shape  was  elliptic.  Dugan  et  al  [ 13]  showed  that  the  region  shape 
was  not  elliptic  after  about  one  half  of  a  Brunt-Vaisala  period  in 
the  fully  mixed  case.  However,  the  energetics  are  integrals  of  the 
solution  so  they  tend  to  be  only  weakly  dependent  upon  the  details 
of  the  mixed  region  shape.  The  validity  of  the  assumption  of  an 
elliptic  shape  was  shown  by  Dugan  et  al  [ 13]  by  computing  the  ener¬ 
getics  numerically  for  times  up  to  about  one  Brunt-Vaisala  period. 
The  results  agreed  reasonably  well  with  those  predicted  by  Bell  and 
Dugan  [11].  The  mixed  region  shapes  predicted  numerically  here  and 
shown  in  Figure  2  are  even  more  elliptic  than  in  the  fully  mixed 
case  so  the  assumption  of  an  elliptic  shape  for  energetics  computa¬ 
tions  is  a  better  one  here. 


From  the  above  discussion,  it  is  apparent  that  the  energy 
content  of  the  mixed  region  can  be  computed  from  its  width.  For  an 
elliptic  shape  in  the  uniform,  partially  mixed  case,  the  potential 
energy  is 


(8a) 


valid  to  the  Boussinesq  approximation,  and  the  kinetic  energy  is 


(8b) 
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The  values  of  a(0  and  d<Vtft  can  be  obtained  from  Figure  3  to  give 
the  results  shown  in  Figure  5.  The  case  of  £=.Qf  is  illustrated. 


1  2  3  4  5  6 

Nt 

Figure  5 

The  result  is  sufficient  to  conclude  that,  as  in  the  fully  mixed 
case,  the  energy  content  of  the  mixed  region  is  largely  given  up  to 
the  exterior  fluid  in  one  Brunt-Vaisala  period.  This  case  of  €=•,©/ 
is  quite  representative  of  the  other  caseg  run  so  plots  of  those 
results  are  omitted. 

Figure  5  also  shows  the  kinetic  and  potential  energy  content 
of  the  mixed  region  as  predicted  by  equations  (8a, 8b)  in  which  ott) 
and  d^/ci'fc  are  given  by  the  theory  of  Hartman  and  Lewis  [7],  i.e. 
expression  (6)  here.  Although  the  result  for  the  kinetic  energy  is 
somewhat  different  and  the  linear  theory  strictly  is  not  valid  in 
this  case,  the  results  predicted  by  the  linear  theory  support  the 
numerical  prediction  of  the  quick  release  of  energy  by  the  mixed 
region . 
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Gaussian  Mixing 


The  density  perturbation  of  Gaussian  form  in  case  II  is  quite 
diffuse  compared  with  the  case  I.  Analyses  of  mixed  region  shapes 
and  widths  are  not  very  meaningful  since  the  edge  of  the  mixed  region 
is  not  clearly  defined.  Nevertheless,  approximations  to  shapes  can  be 
determined  if  one  is  content  to  follow  fluid  particles  that  have  some 
initially  given  positions. 

The  shape  of  a  dense  mass  of  particles  that  is  initially  in 
a  circle  having  the  radius  00  is  plotted  with  time  in  Figure  7. 
Compared  with  the  previous  case,  there  is  an  interesting  difference 


TIME  0 


1.0 


TIME  IN  MRU  NT -VAISALA  PERIODS 


Figure  7 

in  the  final  equilibrium  positions  of  the  particles.  The  shape  is  not 
elliptic  at  all  but  it  shows  a  distinct  peak  near  the  original  vertical 
centerline.  This  equilibrium  shape  is  explained  by  the  type  of  density 
perturbation.  Since  the  amount  of  mixing  goes  to  zero  smoothly  near 
the  edge  of  the  mixed  region,  the  final  equilibrium  positions  of  the 
particles  near  the  top  and  the  bottom  of  the  region  are  not  very  far 
from  the  initial  positions.  Thus,  although  the  particles  near  the 
center  move  appreciably  because  of  the  greater  amount  of  mixing  there, 
those  near  the  upper  and  lower  edge  move  only  slightly. 

The  energy  content  of  the  mixed  region  is  not  clearly  defined 
since  there  is  no  'edge'  to  the  mixed  region  so  no  computations  are 
made . 


Modified  Gaussian  Mixing 

As  in  the  previous  case,  numerical  results  have  been  obtained 
for  the  case  III.  The  particle  distribution  that  was  initially  in  the 
region  Ok  is  displayed  versus  time  in  Figure  8.  This  initial  density 
perturbation  is  closer  to  the  case  I  having  a  sharp  edge  but  the  small 
amount  of  mixing  near  the  edge  in  this  case  continues  to  give  final 
equilibrium  shapes  with  a  peak  near  the  vertical  centerline. 
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Finally,  one  more  observation  about  the  final  width  of  the 
mixed  region  is  useful.  The  potential  energy  for  a  mixed  region  of 
elliptic  shape  is  given  by  expression  (8a).  Since,  in  the  final 
stage  of  collapse,  the  potential  energy  is  all  given  up  to  the 
exterior  fluid,  this  expression  implies  that  the  final  equilibrium 
half-width  of  the  mixed  region  is 

G.  =  CXoe  1  (9) 

This  expression  is  in  good  agreement  with  the  numerical  results  as 
shown  in  Figure  6. 


€ 


Figure  6 

The  result  is  compared  with  that  of  Hartman  and  Lewis  [7]  which 
implies  that 


O  c  Q0  C 


(10) 


and  which,  evidently,  is  valid  only  in  the  limit  of  very  small  mixing. 
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The  final  equilibrium  shape  of  partially  mixed  regions  is  shown 
to  depend  crucially  upon  the  initial  density  perturbation  that  defines 
the  mixed  region.  The  equilibrium  shapes  obtained  for  uniformly  mixed 
regions  are  elliptic  in  cross  section  but  the  shapes  obtained  in  cases 
of  nonuniformly  mixed  regions  show  a  distinct  peak  near  the  vertical 
centerline . 


The  case  of  uniform,  partial  mixing  is  less  realistic  than  the 
nonuniformly  mixed  cases  but  it  lends  itself  to  energetics  computations. 
In  this  case,  for  any  degree  of  mixing,  the  energy  content  of  the  mixed 
region  is  practically  all  given  up  to  the  exterior  fluid  by  the  time  of 
one  Brunt-Vaisala  period  after  the  initiation  of  collapse.  There  is 
little  reason  to  believe  that  this  result  would  be  any  different  for 
other  types  of  mixing.  Also,  the  final  equilibrium  width  of  the  mixed 
region  is  close  to  the  width  predicted  by  a  simple  theory;  that  is 

-<  — *  oo  °  °  v  J 

where  k\0  is  the  Brunt-Vaisala  frequency  outside  the  region  and  N-  is 
the  initial  uniform  Brunt-Vaisala  frequency  inside  the  region. 


The  temporal  dependence  of  the  mixed  region  width  in  the  case 
of  uniform,  partial  mixing  is  compared  with  analytic  results  from  an 
integral  energetics  model  and  a  linearized  differential  model.  Non¬ 
linear  effects  are  shown  to  be  important  in  mixed  regions  that  are 
mixed  as  much  as  207.  or  more;  that  is,  when  ^  . 
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The  purpose  of  this  section  is  to  give  details  on  the 
additional  numerical  features  introduced  by  the  stretching  of  the 
mesh  employed  in  Dugan  et  al  [l3~  in  both  the  x-  and  ?-  directions. 
The  reader  is  referred  to  section  2  of  the  above  reference  to  obtain 
a  complete  description  of  the  conservation,  accuracy  and  stability 
of  the  constant  mesh  scheme;  a  discussion  is  also  contained  in 
Williams  '15  and  Piacsek  and  Williams  f  16\ 

Figure  9  illustrates  the  stretched  mesh  in  one  direction  and 
the  arrangement  of  the  variables  on  their  respective  grid  points. 

A  discussion  of  the  general  approach  to  conservation  on  stretched 


Figure  9 

meshes  is  given  by  Bryan  ' 17,  and  will  not  be  repeated  here.  However 
we  will  illustrate  below  the  essential  ideas  by  applying  the  sctieme 
to  the  one-dimensional  transport  of  momentum  and  temperature.  For 
inviscid  flows,  the  equations 

H  “  ~  a**  "  *x  (1.1) 

vr  _  vr 

'~a^  (1.2) 

will  serve  adequately  to  demonstrate  the  technique.  Using  Figure  9 
and  the  quadratica 1 ly  conserving  scheme  of  Piacsek  and  Williams  [16], 
we  may  write 


Note  from  Figure  9  that  and  &*•  refer  to  elements  of  different 

mesh  sequences,  one  'centeredron  and  the  other  onT  , 

respectively.  Upon  multiplying  expression  (1.3)  by  AX  .  y  and 
expression  (1.4)  by  and  summing  over  i,  the  sums  encountered  in 
Dugan  et  al  13  are  produced  and  these  were  shown  to  be  conserved 
However,  whereas  the  T-boxes  are  actually  centered  on  Tt  ,  the  u-boxes 
are  centered  on  U  ivl<v ,  resulting  in  a  loss  of  second  order  accuracy 
in  the  truncation  error  unless  the  mesh  is  stretched  in  certain 
special  ways.  The  current  mesh  stretching  has  only  first  order 
accuracy,  and  is  given  by 

,  V,c  v,  *Q.  -L 

^  (  e  +1)  (1.5) 

In  the  calculations  in  this  repor*-,  the  mesh  was  stretched  such  that 
t.ie  final  grid  length  in  the  horizontal  was  eighteen  times  the 
initial  length,  and  the  total  box  size  was  eight  radii  wide  and 
four  radii  high.  The  stretching  enabled  a  reduction  in  the  number 
ot  horizontal  mesh  points  by  a  factor  of  five. 

In  addition  to  the  changes  outlined  above,  the  computations 
performed  in  solving  the  Poisson  equation  in  Dugan  et  al  [ 13]  are 
now  done  in  double  precision  so  that  roundoff  errors  are  minimized. 

is  gives  an  order  of  magnitude  improvement  in  the  conservation  of 
total  energy  in  the  computation  box;  the  change  in  total  energy  in 
fourteen  Brunt-Vaisala  periods  of  computation  is  less  than  .57. 
versus  somewhat  less  than  57.  in  single  precision  arithmetic. 


